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We have applied the Fast Fourier transform (FFT), which allows to compute efficiently convolution 
sums, to solve the set of self-consistent T-matrix equations to get the Green function of the two 
dimensional attractive-U Hubbard model below T c , extending previous calculations of the same 
authors. Using a constant order parameter A(T), we calculated T c as a function of electron density 
and interaction strength U. These global results deviate from the BCS behavior remarkably. 
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Although the Hubbard modeE is the most simple model to describe correlated electron behavior in a solid, the 
mathematical treatment is far from trivial. Many attempts have been made to understand the phase diagram. A 
fully understood Hubbard model might form the basis of understanding correlated electron systems as much as the 
Ising model did for understanding critical phenomena in magnetism. The attractive-U Hubbard model might play 
an important role in understanding high-temperature superconductivity and has been attracting much attention in 
the past few years. We have implemented the T-matrix approximation which goes beyond the usuaUmean-field 
approximation and becomes exact in the dilute limit, i.e. where only two-particle interactions take placai- 

We consider the attractive-U Hubbard model in two dimensions on a square lattice (lattice constant ap 

H = J2^ C L C ^ + U H C k+qT C -kl C -k'l C k' + qT (1) 
k.cr k,k',q 

with band energy ek = —2t(cos h x &-\-coshy(i) and on site attraction U <C 0, £k — ^k M) M is the chemical potential, and 
t, the hopping of electrons between nearest neighbour sites, determines the energy unit. The creation (annihilation) 
operators for an electron with momentum k and spin a are denoted by cj^ (cko-)- 

The T-matrix (effective interaction) is the sum of particle-particle ladder diagrams with the smallest number of 
closed fermion loops. In the low-density limit, these are the dominating terms of the perturbation expansion in terms 
of the interaction U. For the Hubbard model, where we have only on-site interactions, the T-matrix approximation 
leads to a set of self-consistent equations for the one-particle Green functionu in the normal phase 

G(k,iu n ) = [iuj n - ^ k + T,(k,iu n )]-\ (2) 

where the diagonal self-energy term 

E(k, iuj n ) = — ^ r (q> «£m)G(q - k, ie m - iu n ), (3) 

q,m 

depends on the T-matrix 

T(<l,ie n )= ~V , (4) 
1 + C/x(q,«e«) 
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which is a simple function of the independent pair-susceptibility 



X(q, ie n ) = G(k, iw m )G(q - k, ie n - iu m ). (5) 



k,m 



Here, uj n = (2n + 1)5 and e„ = 2nS are the fermionic and bosonic Matsubara frequencies, N = N x N y where N x and 
N y are the grid dimensions in k-space and (3 = 1 /T is the inverse temperature. We fix the chemical potential by the 
electron density (one spin direction): 

n(J3, m) = Urn V G(k, ^ n )e iw »". (6) 

k,n 

To go below T c , we introduce a constant order parameter A = |A| into the Green function following the usual 2x2 
Nambu matrix formalism and get for the diagonal part the expression (* means complex conjugate) 

r>(\ ■ \ r*(k,^») 

G(k ^ )= |r(k^ n) p + |Ap (7) 

where 



T(k, iu n ) = iiv n - £k + E(k, iu) n ). (8) 

Equation (Q) reduces to the usual BCS Green function when the self-energy term is set to zero (or to the Hartree 
shift). A is our approximation for Si2(k, iw n ), the off-diagonal self-energy. The order parameter A(T) is determined 
by 

I = J-V I (9) 

U (3N^ |r(k,i Wn )l 2 + |A| 2 {> 

which closes the set of equations. 

To solve the set of equations (||), (Q), (^), (^), (R) and (||) we apply the following scheme: 

1. Start by calculating Go(k, iu) n ), i.e. the Green function for the free system (E = 0). A suitable initial value for 
/j, and A must be given (/i = —3.5 and A = 0.5 are reasonable values for T = 0.1, U = —4 and n = 0.1). 

2. Calculate xi^^m), T(q,ie m ) and S(k, ioj n ) using equations (||), ffl) and (||). 

3. At this point, we need an improved estimation of /i and A. To get a stable iteration scheme, both parameters 
must be adjusted simultaneously. We do that by searching a solution (//, A) for equations (||) and (||) using a 
Newton algorithm. For technical purposes we neglect the dependence on fj, and A for the self-energy to be able 
to numerica 

lly calculate the partial derivatives needed for the Newton algorithm. Because of that approximation we can not 
take the new estimate of (//, A) directly. Instead we go only a small step, in the fj, — A plane, from the current 
point to the new point (about one third of the total distance) . 

4. Calculate an improved Green function using the new parameters /i and A. 

5. Repeate steps |]to |] until the electron density has reached its desired value within a given tolerance. 

In order to obtain results which are independent of finite size, one should use at least some 10 3 Matsubara frequencies 
and a grid of 30 x 30 lattice points. The above scheme works in principle but a closer look at the equations for \ 
and £ shows that the straight forward implementation of these equations does not work in practice. This is due to 
the 4-fold loops which would occur in the computer program. Suppose we use 2000 Matsubara frequencies and a 
30 x 30 grid. Then we have to carry out for every grid point and every frequency the double sum over all frequencies 
and all grid points. This are of the order (30 2 x 2000) 2 = 3.24 x 10 12 complex operations. Even one of the fastest 
super-computers would need one to several hours to make one iteration step. 

Since the frequency and momentum summations are convolutions, we evaluate them using the fast Fourier transform. 
The transforms k — * x and k <— x are the usual ones and we do not elaborate them any further. The transforms from 
t — > iuj and t <— iu> arc described in more detail. In the following, the notation 
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M-l 

M . 
and 

M-l 

'M 



M-l 

TTT M [F( Xj )] n = ^=Y1 e- 2 ^F( Xj ) (10) 



M-l 

^T- 1 [F(x„)] j = -= £ e 2 ^F(x n ) (11) 



n=0 



is used. The Matsubara frequencies are slightly redefined to be more suitable (having non-zero indices) for numerical 
work and read: 



We discretize the integral 



w„= (2n + l-M)^, e n = (2n-M)^. (12) 



r/3 

= / dre^" T G(r) (13) 



by writing tj = j At, j = . . . M — 1 where At = and M is the number of Matsubara frequencies used. We obtain 

G(^„) = ^OT-y^-^Gfo)],, (14) 



The phase factor e 1 "- 7 ^ m arises because of the fermionic frequencies u> n and the shift in the definition of oj„ . Similarly 
one can define the transforms for the bosonic frequencies as 

X(ie n ) = -j=TTTm [e^'Jf (r 3 -)]n, (15) 
vM 



where X is either \ or S. We can now rewrite equations (||) and (|^) to read: 



X (k„) = -^=OT^[e-^G 2 (r J -)]„ (16) 
v A/ 



and 

£(z„„) = -^OT M 1 [ e ^^- 1 )T(T,)G(-T,)]„. (17) 

These expressions are also well suited for parallel machines since the 3D-FFT (two space and one imaginary time 
dimension) can be decomposed into parallel processes. 

We also need to calculate the electron density, but evaluating (||) numericaly is not possible. To remove the limit 
77 — s- + , we make use of the definition of the Green function 

G(t) = (ct(r)c(O)), t>0 (18) 

and 

G(-t) = (c(0)c t (-T)) = -1 + (J(-t)c(0)), t>0. (19) 
Taking the sum G(r) + G(— r) and letting r — > + we get 

„=i[G(0+) + G(0-)] + i = G(0) + i. (20) 



Therefore, the electron density now reads 

'R it \ = 

2 (3N 
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an expression which can be evaluated easily. A similar correction for G(t — 0) must be applied in equations ( |l6| ) and 
(U). Our Eq.@ generalizes Eq.(3.1.2) of Mahan'sH. 

We have compared n((3,fi) for the cases E = and E = nU/2 (the Hartree shift) with the exact results and 
concluded that a grid of 32 x 32 lattice sites and 2048 frequencies is a lower limit for U = —4 and temperatures down 
to T = 0.1. For larger \U\ and/or lower temperatures one should increase the number of frequencies. 

So far we have evaluated the Green function at the Matsubara points. Normally, we are really interested in the 
Fourier transform of the retarded real-time Green function which is a function of real frequency. In principle, we 
get this function by analytic continuation of the complex frequency Green function from the Matsubara points to 
the real axis. Since the resulting integral equations would be more complicated (involving integrals over Fermi and 
Bose distribution functions) than the discrete frequency summations, we first calculate the Green function at the 
Matsubara points as described in the preceeding section. Then we coatinue to the real frequency axis by fitting a 
rational function (M-point Pade approximant) to the calculated valucsa. The dynamical properties of the attractive 
Hubbard model in the superconducting phase is discussed in Ref.Q. Here, we concentrate on the global properties of 
the attractive Hubbard model. 

The algorithm works as follows: Given a function f(zi) = Ui with values u% at M complex points z%, i = 1,2,..., M, 
the Pade approximant is defined as a ratio of two polynomials which can be written as a continued fraction 

C »W = a^-zQ ™ 
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a M {z - z M -i) 



where the coefficients a* are to be determined so that 

C M (zi)=Ui, i = l,2,...,M (23) 
which is fulfilled when the aj are given by the recursion 

ai=gi(zi), gi(zi) = Ui, i = 1, 2, ...,M (24) 

and 

[z - z p -i)g p -i{z) 

Once the coefficients ai are determined for a particular function, the function values at a real frequency ui can be 
obtained by setting z = u + iS in (|22|), where S is used to remove unphysical peak structures due to finite size effects. 
The choice 6 = 0.1 . . . 0.2 is appropriate for our calculations. In reality, the choice of S is given by the precision avith 
which the self-consistency in the Matsubara frequencies is done. This point has been discussed by Georges et aH. In 
this paper we limit the discussion of the results to the regime T < T c , since results for T > T c are reported elsewhereS. 

Figure |l| shows the temperature dependence of the order parameter for n = 0.1 and U = —At. The system size is 32 
by 32 grid points in k-space and 1024 Matsubara frequencies. The order parameter shows a very sharp drop to zero 
(compared with BCS behavior) which indicates the strong influence of fluctuations for all temperatures. A(T = 0) is 
nearly twice as large as the corresponding A bcs whereas the critical temperature is much lower then T^ cs by more 
than a factor three. For example, A(0)/T c as 3.86, which is more than twice the BCS ratio (si-,1.76). This implies 
that we are not in the weak coupling limit, opposite to the case of Martin- Rodero and FloresOJJ who find the BCS 
universal ratio in second order of perturbation for the continuous Hubbard model. The value of A(T) = defines 
the critical temperature, T c . In order to calculate T c we have drawn a straight line close to the sharp drop of the 
order parameter. Near the critical temperature we had some convergence problems. However, the evaluation of T c is 
precise because it is calculated as described previously. 

The density dependence of T c is plotted in Figure ^. It is well known, that in the strong coupling limit the 
attractive-U Hubbard model close to half-filling can be mapped to a Heisenberg model which shows no Kosterlitz- 
Thouless transitiorU. It is therefore expected, that the critical temperature near half-filling is reduced even in the 
weak- and intermediate-coupling regime. The lack of this property in our results can be attributed to the neglect of the 
particle-hole channel (charge fluctuations) in the T-matrix. Then, we should restrict ourselves to low densities where 
the T-matrix approach is indeed valid. To treat higher densities, we should implement, for example, the FLEXC 
approacht3 goal which we are pursuing at this moment for the attractive Hubbard model. Here, w£_mention that 
the T-matrix and FLEXC approaches are conserving in the Kadanoff-Baym sensso. Denteneer et alE£l obtained the 
critical temperature by calculating the helicity modulus associated with a wavelike distortion (Aj = |A|exp(2iqr ■)) 
of the order parameter in the BCS approximation as a function of temperature and subsequent comparison with the 
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Kosterlitz-Thouless relation between critical temperature and helicity modulus. But they miss the logarithmic drop 
to zero when approaching half-filling. 

Figure || shows the critical temperature as a function of interaction strength. In contrast to the BCS behavior 
(which shows an exponential increase of T c small \U\ and goes linear in \U\ for larger values), the increase of T c is 
reduced drastically. The expected decrease of T c for large \U\ can not be observed for the values of U treated here. 
The expectation is based on the fact, that for strong attraction and near to half-filling the model can be mapped 
onto a pseudo-spin model with effective interaction constant J = 4t 2 /\U\, resulting in a T c decreasing with increasing 
\U\. The behavior we observe for T c /t vs U/t is similar to the one obtained by Nozieres and Schmitt-Rinlo using 
the Thouleas' criteria in the normal phase. Most likely we have to include additional fluctuations in the off-diagonal 
self-energylij. Contrary to the normal state calculations where, for p — 0.1, we could not obtain convergence for 
\U\ > 4.0, now in the superconducting phase the program is stable for large values of \U\. 

In conclusion, we have calculated the critical temperature of the negative-U Hubbard model within the T— matrix 
approximation. The expected smooth cross-over from BCS to Bose-condensation can not be fully observed in the 
parameter space studied in this paper. Although the exponential-to-linear increase of T® cs with growing \U\ is 
reduced drastically, an optimal critical temperature has not been found. PerhapSr-we should fullyj-go beyond BCS in 
the off-diagonal self-energy, as it has been done by Shafroth and Rodriguez-Nufieai3. In this worktLa, the authors have 
studied the dynamical properties of the attractive Hubbard model in presence of double fluctuations. As previously 
discussed, we have obtained that A/T c is more than twice the BCSLyalue. To get T c = at half-filling we must include 
charge fluctuations. Investigations along these lines is under waj£3. The value A(0)/T c w 3.86 and the temperature 
behavior of A(T) at U/t = —4.0, n = 0.1 are very different with respect to BCS results. These global results suggest 
that for U/t — —4.0, n = 0.1 we are already in the intermediate coupling regine where correlations are important. 
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FIG. 1. Order parameter A versus temperature for density n = 0.1 and interaction strength U = —4. 
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FIG. 2. Critical temperature T c versus electron density n (n = 0.5 corresponds to half-filling) for an interaction strength 
U = -4. 
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FIG. 3. Critical temperature T c for various interaction strengths U for constant density n = 0.1. 
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